#include <iostream>
#include <vector>
#include "Geometry.h"
#include "Excitation.h"
#include "Solution.h"

int main()
{
	clock_t START,END;
	START = clock();

	std::vector<TriangleMesh*> stuff;
	const int N = 1;
	//Import PO region
	Object* PO = new Object();
	PO->ReadGeometry("5m_plate.dat",emx::Real3D(0,0,0));
	stuff.push_back(PO);

	//Import Originals
	Geometry* Origin = new Object();
	Origin->ReadGeometry("cube.dat",emx::Real3D(0,0,0));

	//Import Equivalences
	EquivalentSurface* Equivalence[N];
	Equivalence[0] = new EquivalentSurface();
// 	Equivalence[1] = new EquivalentSurface();
// 	Equivalence[2] = new EquivalentSurface();
// 	Equivalence[3] = new EquivalentSurface();
// 	Equivalence[4] = new EquivalentSurface();
// 	Equivalence[5] = new EquivalentSurface();
// 	Equivalence[6] = new EquivalentSurface();
// 	Equivalence[7] = new EquivalentSurface();
	Equivalence[0]->ReadGeometry("sphere_mid_sparse.dat",emx::Real3D(0,0,0));
// 	Equivalence[1]->ReadGeometry("sphere_mid_sparse.dat",emx::Real3D(-2,2,-2));
// 	Equivalence[2]->ReadGeometry("sphere_mid_sparse.dat",emx::Real3D(2,-2,-2));
// 	Equivalence[3]->ReadGeometry("sphere_mid_sparse.dat",emx::Real3D(2,2,-2));
// 	Equivalence[4]->ReadGeometry("sphere_mid_sparse.dat",emx::Real3D(-2,-2,2));
// 	Equivalence[5]->ReadGeometry("sphere_mid_sparse.dat",emx::Real3D(-2,2,2));
// 	Equivalence[6]->ReadGeometry("sphere_mid_sparse.dat",emx::Real3D(2,-2,2));
// 	Equivalence[7]->ReadGeometry("sphere_mid_sparse.dat",emx::Real3D(2,2,2));
	stuff.push_back(Equivalence[0]);

	//Set Excitation
	Excitation* Ex = new PlaneWave(emx::C, 0, 0 ,0);
	Solution* proj = new MoM_Solution(Ex);
	
	//Save the Incident Wave
	std::vector<emx::zVector> Inc;
	Inc.resize(stuff.size());
	for (int i=0; i<stuff.size(); i++){
		Inc[i].resize((stuff[i]->GetMedium() == emx::PEC)?stuff[i]->GetEdgeNum():2*stuff[i]->GetEdgeNum());
		proj->Fill_inc(stuff[i],Inc[i]);
	}

	//Calculating Matrix Elements
	emx::zMatrix S(2*Equivalence[0]->GetEdgeNum(),2*Equivalence[0]->GetEdgeNum());
	proj->ExtractMatrixS(Equivalence[0],Origin,S);

	blitz::Array<emx::zMatrix,2> T;
	T.resize(stuff.size(),stuff.size());
	clock_t start,end;
	std::cout<<"Filling Matrix T..."<<std::endl;
	start = clock();
	omp_set_num_threads(omp_get_num_procs());
#pragma omp parallel for
	for (int i=0; i<stuff.size(); i++){
		for (int j=0; j<stuff.size(); j++){
			if (i != j){
				int size1,size2;
				if (emx::Vaccum == stuff[i]->GetMedium()){
					size1 = 2*stuff[i]->GetEdgeNum();
				}
				else{
					size1 = stuff[i]->GetEdgeNum();
				}
				if (emx::Vaccum == stuff[j]->GetMedium()){
					size2 = 2*stuff[j]->GetEdgeNum();
				}
				else{
					size2 = stuff[j]->GetEdgeNum();
				}
				T(i,j).resize(size1,size2);
				proj->ExtractMatrixT(stuff[i],stuff[j],T(i,j));
			}
		}
	}
	end = clock();
	std::cout<<"	Filling Matrix T consumed "<<double(end-start)/CLOCKS_PER_SEC<<" seconds."<<std::endl;

	std::cout<<"Iterating..."<<std::endl;
	start = clock();

	//Solve the Problem Recursively
	for (int k=0; k<3; k++){
		for (int i=0; i<stuff.size(); i++){
			proj->Inc_Initial(stuff[i],Inc[i]);
			for (int j=0; j<stuff.size(); j++){
				if (j != i){
					proj->Transfer(stuff[i],stuff[j],T(i,j));
				}
			}
			if (stuff[i]->GetMedium() == emx::Vaccum){
				proj->Fill_sca(stuff[i],S);
			}
		}
	}
	end = clock();
	std::cout<<"	Iterating consumed "<<double(end-start)/CLOCKS_PER_SEC<<" seconds."<<std::endl;

	//Calculating the RCS
	proj->GetResults(&stuff[0], stuff.size(), 0,180,181,0,0,1,0);

	END = clock();
	std::cout<<"\n\nThe complete solution consumed "<<double(END-START)/CLOCKS_PER_SEC<<" seconds."<<std::endl;
	getchar();

	return 0;
}